Role of fourth-order phase-space moments in collective modes of trapped Fermi gases 
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We study the transition from hydrodynamic to coUisionless behavior in collective modes of ul- 
tracold trapped Fermi gases. To that end, we solve the Boltzmann equation for the trapped Fermi 
gas via the moments method. We showed previously that it is necessary to go beyond second-order 
moments if one wants to reproduce the results of a numerical solution of the Boltzmann equation. 
Here, we will give the detailed description of the method including fourth-order moments. We apply 
this method to the case of realistic parameters, and compare the results for the radial quadrupole 
and scissors modes at unitarity to experimental data obtained by the Innsbruck group. It turns out 
that the inclusion of fourth-order moments clearly improves the agreement with the experimental 
data. In particular, the fourth-order moments reduce the effect of collisions and therefore partially 
compensate the effect of the enhanced in-medium cross section at low temperatures. 

PACS numbers: 67.85. Lm,03.75.Ss 



I. INTRODUCTION 

The study of collective modes of trapped two- 
component Fermi gases revealed interesting information 
about different dynamical regimes [l|. Initially, the aim 
was to find signals for the superfluid-normal phase tran- 
sition. However, near a Feshbach resonance, the atom- 
atom scattering cross section can be large enough to en- 
sure (normal-fluid) hydrodynamic behavior of the gas 
above the superfluid critical temperature Tc. In this 
case, a change in the behavior of the gas is observed at 
much higher temperature, when the gas gets more and 
more dilute until the coUisionless regime is reached. The 
most interesting modes in this context are those which 
in the coUisionless case exhibit deformations in both co- 
ordinate and momentum space. Such modes are, e.g., 
the quadrupole and the scissors modes. In the presence 
of superfluidity or collisions, the deformation of the mo- 
mentum sphere is suppressed, so that the frequencies of 
these modes are different from those in the coUisionless 
case. In the intermediate regime, the damping of these 
modes is very strong. Both the radial quadrupole mode 
and the scissors mode were experimentally studied by the 
Innsbruck group [i|-Q . 

From the theoretical side, the continuous transition 
from coUisionally hydrodynamic to coUisionless behav- 
ior can be studied by using the semiclassical Boltzmann 
equation. At present, there is no technique which would 
allow for a fully quantum mechanical description of col- 
lective modes of systems containing a few hundred thou- 
sand particles, including the coUisional effects. But even 
the solution of the Boltzmann equation is far from being 
simple, and most of the time further approximations are 



made. A very common approximation is the relaxation- 
time approximation, which was used, together with the 
so-called scaling ansatz, to describe collective modes and 
the expansion of the gas after the trap is switched off [3l ■ 
In the case of collective modes, this method is equiva- 
lent to the method of phase-space moments up to second 
order, which was applied to the radial quadrupole, scis- 
sors, and breathing modes @, [3 Q- In both methods, 
the phase-space distribution function is constrained to a 
simple form, but the advantage is that one can perform 
computations almost analytically. 

There are other approaches like the fully numerical so- 
lution of the Boltzmann equation as developed, e.g., in 
Refs. [7H10|. In this case no constraint is put on the func- 
tional form of the distribution function, but the price to 
pay is that the computations are very time consuming. 
Maybe the computation time could be significantly re- 



duced by using new adaptive algorithms |ll| . but to our 
knowledge, no numerical calculation has been performed 
so far for degenerate Fermi gases with parameters (num- 
ber of atoms, trap geometry) corresponding to real ex- 
periments. 

In our previous work [9[ , we compared the results of a 
numerical simulation of the quadrupole mode in a spher- 
ical trap containing a reduced number of atoms with 
the corresponding results of the second-order moments 
method. The surprising outcome was that the second- 
order moments method strongly overestimates the effects 
of collisions. This problem could be cured to a large ex- 
tent by generalizing the method of moments to fourth 
order. This is our main motivation for the present work, 
where we apply the fourth-order moments method to the 
radial quadrupole and scissors modes in a realistic trap 



geometry, and compare it directly with the experimental 
results of Ref . Q • 

The Boltzmann equation requires as microscopic in- 
put the mean-field potential (in the Vlasov part) and the 
cross section (in the collision integral). Here, the mean 
field will be neglected since we found in Ref. Q that 
it affects only very weakly the frequencies and damping 
rates of the collective modes near unitarity. The main 
effect comes clearly from the collisions. In the case of 
large scattering length a and temperature slightly above 
Tc, one expects the Fermi gas to be in the "pseudogap 
regime" in which pair correlations play an important role 
although the pairs are not condensed. In this regime, 
the relaxation time is strongly reduced since the scatter- 
ing cross section calculated in the surrounding medium 
is enhanced as compared with the free one [12l | - an ef- 
fect which in the context of nuclear physics has already 
been known for a couple of years |13| . Previous stud- 
ies 0, @ using the in-medium cross section in trapped 
Fermi gases found that this reduction of the relaxation 
time badly deteriorates the agreement with the experi- 
mental results. Here we argue that this discrepancy was, 
at least to some extent, due to the failure of the second- 
order moments method and not due to the enhancement 
of the in-medium cross section. 

The paper is organized as follows. In Sec. [II] we de- 
scribe the method, starting with a very general formu- 
lation and specializing then to the scissors and radial 
quadrupole modes. We explain how the response of the 
system is obtained and how we extract from it the fre- 
quencies and damping rates. In Sec. IIIII we discuss our 
results. We show how the inclusion of fourth-order mo- 
ments affects the response function and the correspond- 
ing frequencies and damping rates and compare the the- 
oretical results with experimental data. In Sec. IIVI we 
summarize and give an outlook to future studies. Some 
technical details are given in the appendix. 

Throughout the article, we use units with ft = fc^ = 1. 



II. FORMALISM 



of collective modes as compared with the effects coming 
from collisions between atoms [6| . 

In the normal fluid phase and under other assumptions 
we already discussed in Ref. [6[ , we can describe the sys- 
tem with a semiclassical distribution function /cr(r, p, i), 
where a =t, |. We restrict ourselves to the case of an 
unpolarized gas (A'^ = N^ = N/2) and to excitations 
where the two components move together: f^ = f^ = f. 
The normalization of / is^ 



(27r)3 



f{r,p,t) = ^ 



(2) 



and the average value of a generic quantity x(r, p) is 
2 /■ d^rd^p ^, ^ , , 

In equilibrium, the distribution function reads 

feqir, P) 



(3) 



g/3[pV2m+V(r)-M] ^ I ' 



(4) 



where /? = l/T is the inverse of the temperature and /i 
is the chemical potential. 

When the system is excited, the time evolution of / is 
governed by the Boltzmann equation IJ]. We consider 
small perturbations Sf of the distribution function from 
equilibrium and write them as 

<5/(r, p, t) = /e,(r, p)[l - /e,(r, p)]$(r, p, t) . (5) 

The function <I>(r,p,t) can be assumed to be smooth 
since the fact that Sf is peaked near the Fermi surface is 
already accounted for by the prefactor feq{^ ~ feq)- The 
linearized Boltzmann equation then reads 



/e,(l - feq) $ + il • V,$ ~ VrV ■ Vp$ 



+ P—-Vr5V]=-I['^]. (6) 

m J 

The linearized collision integral in the right-hand side is 



A. Moments method for the Boltzmann equation 

We consider a two-component (t, ^) gas of Fermi atoms 
of mass m and with an interspecies attractive interaction 
(the scattering length is a < 0). The gas is loaded in a 
harmonic, usually anisotropic, trap 



V{v) 



'/ 2 2 , 2 2 I 2 2\ 

-(w^x -fw^y +u;^z ) 



(1) 



Moreover, since we are interested in collective modes and 
their time evolution, we include in the external potential 
felt by the atoms a small time-dependent part (5F(r,i), 
that will be used to simulate the excitation of the mode. 
As mentioned before, the mean field felt by the atoms 
due to their interaction will be neglected here, since at 
unitarity it is only of minor importance for the properties 



/[$] 



(27r)3 



d?L m ■'"i-'^'i^ 



X (1 - /' )(1 - /' i)($ + $1 - *' - $;) . (7) 



The various feq and $ are all evaluated at the same r, t 
but at different momenta p, p^, p', or p[, respectively, 
which due to momentum and energy conservation satisfy 
P + Pi = P' + P'l = k and |p - Pi| = |p' - pi| = 2g. 
The solid angle between the initial and final relative 
momenta in the center-of-mass frame, q and q', is de- 
noted O. The cross-section da/dfl used in the present 



^ Notice that this normalization diflfers from that given in Ref. [a| 
by a factor {2tt)^. 



paper is the in-medium cross section which is calcu- 
lated as described in Ref. [y|. At temperatures close 
to the superfluid transition temperature Tc, this cross- 
section is strongly enhanced with respect to the free one 
dao/dfl = a^/[l + (qa)'^], at least for collision partners 
near the Fermi surface with zero total momentum. 

Since the function $ is supposed to be smooth, one can 
try to approximate it by a polynomial in the components 
of r and p with time-dependent coefficients Ci, 



$(r,p,i) = ^c,(i)0,(r,p) 



(8) 



The choice of the basis functions (j>i depends on the mode 
one wants to describe (sec discussions in Refs. (g. [l5|). 
However, let us first explain the general idea before focus- 
ing on the examples of the radial quadrupole and scissors 
modes. 

In order to obtain the so-called response function, it is 
sufficient to consider a perturbation which is a (5 pulse. 



I.e. 



SVir,t) ^6{t)V{r). 



(9) 



Then the Fourier transform of Eq. ^ with respect to t 
gives 



n 
^Ci(w) feq{l-feq)[ 
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pH^i 
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-/e,(l - /e,)/3- • VrV{v) , (10) 



where c^ {u) is the Fourier transform of c^ (t) . Now we take 
the moments of Eq. ([TU]), i.e., we multiply it by each of 
the basis functions 0i and integrate over phase space. In 
this way, we obtain n coupled linear algebraic equations 
for the n coefficients Ciiuj). In matrix form, they can be 
written as 






where 



Aij = -iujMij 
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and 
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(11) 

(12) 
(13) 



A. 



^./e,(l~./e,)0.{0„|^ + ^}, (14) 



{2^f 



-0i/[0j] , 



(15) 



Oi = 



^^0,/,,(l-/,,)p.V.nr). (16) 



The contribution ^*™"* of the transport part of the 
Boltzmann equation to Aij has been written in a com- 
pact form using the Poisson brackets {•,•}. One can show 
that M and ^'=°" are symmetric matrices, while ^*™"-' is 
antisymmetric. 

Once we have solved Eq. ([TT|) for the coefficients Ci(a;), 
we know the time-dependent distribution function fpq + 
6f and we can obtain the time evolution of the average 
of any dynamical quantity using Eq. ^. 

In summary, making a polynomial ansatz for the time- 
dependent distribution function, we reduced the lin- 
earized Boltzmann equation from an integro-diffcrential 
equation to a system of n coupled linear algebraic equa- 
tions for the coefficients c;. 



B. Scissors and quadrupole modes 

Consider an elongated trap with elliptic transverse sec- 
tion (i.e., Wj. > ujy 3> W2) containing a gas in equilibrium. 
The scissors mode is a collective mode that is excited by 
tilting the trap by a small angle (~ 5°) around the z- 
axis. After this excitation, the cloud is rotating back 
and forth around the z axis, and what is measured is the 
time dependence of the angle of the orientation of the os- 
cillating cloud with respect to the trap potential. For the 
details on the experimental realization of this mode and 
the results at finite temperature and different scattering 
lengths, sec Refs. [l|, Q. 

If the initial potential is harmonic, the scissors mode 
is excited by the perturbation 



V{r) — axy . 



(17) 



where a is a factor characterizing the strength of the per- 
turbation. Under the assumption that the shape of the 
cloud does not change during the oscillation, the mea- 
sured angle is proportional to the expectation value 



Q{t) = (xy) 



(18) 



The minimal ansatz for the function $ that can re- 
Droduce the scissors mode contains four terms and reads 
lii 



$2«d = ClXy + C2PxPy + csxpy + c^ypx 



(19) 



All four terms are of second order in the components of 
r and p. In fact, at second order, there are no other 
combinations which satisfy the symmetry of this exci- 
tation which is odd under {x,Px) — > {—x,—px), odd 
under {y,Py) — > {—y,—py), and even under {z,pz) — ?> 
(—z,—pz). In the present case of a harmonic potential 
without mean field, this set of basis functions is closed 
with respect to the operators that arc in the transport 
part of the Boltzmann equation, i.e., on the left-hand 
side of Eq. ©. 

As noted in Ref. [0| in the case of the quadrupole mode 
in a spherical trap, the method of second-order moments 



strongly overestimates the collisional effects because it 
implicitly neglects the position dependence of the relax- 
ation time r. Remember that the effect of collisions is to 
produce hydrodynamic behavior by maintaining the mo- 
mentum distribution spherical. The deformation of the 
momentum distribution is described by the second term 
in $2rad, i-e., the term ex PxPy So, the corresponding co- 
efficient C2 is large in the case of few or no collisions and 
small in the case of frequent collisions. In the trapped 
system, however, the collision rate is very different de- 
pending on the position: Near the center, the density 
and thus the collision rate is much higher than at the 
surface. Therefore, the deformation of the momentum 
distribution should depend on the position. This can- 
not be accomplished with the ansatz (fTQ]). since the term 
ex PxPy is independent of r. 

Let us therefore go to the next higher order, which is 
fourth order. At this order, terms like x^PxPy etc. appear 
which allow us to describe the position dependence of the 
deformation of the momentum distribution. Keeping all 
terms which respect the symmetries mentioned above, we 
must then include 32 terms into the ansatz for $: 



32 



$4t?i(r,P,0 = ^c^{t)(l),{v,Tp) 



(20) 



The basis functions (pi can be compactly defined in the 
following way: 

0j+4O-i)(r,P) =.g«(r,P)/ij(r,P), (21) 

where i = 1, . . . , 4 and j = 1, . . . , 8, and 



9i=xy, 
hb^pl: 



92 = PxPy 
h2 ~ X . 



he, 



Pv 



, 53 = XPy , .94 = yPx 

hr ^pI, hs^ zp^ . (22) 



It is easily seen that the first four terms of ^4th reproduce 
^2nd, while the subsequent ones are fourth-order terms 
in the components of r and p. 

Let us now turn to another mode, the radial 
quadrupole mode in an axially symmetric trap, lUx = '^y 
In Refs. [3|, Ig, the corresponding perturbation was writ- 
ten as V (X x"^ — y"^ and the measured observable was 
(x^ — y^) . However, since the trap is axially symmet- 
ric, we can rotate the coordinate system by 45° around 
the z axis without changing anything. By doing so, 
one sees immediately that the perturbation is then of 
the form V oz xy and the measured observable becomes 
(xy), like for the scissors mode. In conclusion, the radial 
quadrupole mode is a special case of the scissors mode in 
the limit of equal trap frequencies uJx = uJy , and it there- 
fore does not require any additional effort to describe 
both modes. 



C. Response function 

As already mentioned, we follow the observable Q = 
(xy), which, with our choice of basis functions, can be 



written as Q = {(jii). Using Eqs. dU and ([8|), this expec- 
tation value can be expressed in terms of the coefficients 



Ci as 



2 ^^ 

4=1 



(23) 



where Mn are the elements of the first row of the matrix 
M defined in Eq. (fTSj) . 

Also the vector a on the right-hand side of the linear 
system of equations (fTTj) for the coefficients Ci{(jj) can 
be expressed with the help of the matrix M . Note that 
p • VrV{v) = a{xpy + ypx) = a{<p3 + ^4), so that Eq. 
(ITBl) becomes 



a/3 



(A-f,3 + M,i) . 



(24) 



Now, the linear system of equations pT|) for the co- 
efficients Ci has to be solved. After some algebra (see 
appendix), the result for the response function can be 
written as 



QH 



-2iaP " {MP)^u[{P-^)m + {P-^)ua] 



fe=i 



w - Wfc -I- iTfe 



, (25) 



where V^ + i^k is the A:th eigenvalue of the matrix 
A/~^(^*™"'' -I- ^<=°'') and P is the matrix containing in 
its columns the corresponding eigenvectors. 

It should be pointed out that it is a very tedious work 
to calculate the elements of the matrices Af , ^*™"-' and 
jicoii corresponding to the fourth-order moments. Here, 
we made use of the Mathematica software to derive the 
expressions. After that, the actual numerical calculations 
are quite fast, the only time-consuming part is the Monte- 
Carlo integration of the moments of the collision term in 
j^coii ^ The numerical inversion and diagonalization of 
a 4 X 4 (second-order method) or 32 x 32 (fourth-order 
method) matrix does not pose any problem. More details 
about the calculation of the matrices are given in the 
appendix. 

For the discussion, the imaginary part of (3(0;) is par- 
ticularly useful, since this so-called strength function 
describes the excitation spectrum corresponding to the 
mode under consideration. 



D. Frequencies and damping rates 

In the previous literature (3|, |6[ , where the second-order 
moment method was used, the frequencies and damping 
rates of the collective modes were identified with the real 
and imaginary parts of the solutions of the characteristic 
equation detA = 0. These are of course equal to the 
imaginary and real part of the eigenvalues of the matrix 
^-i(^trans_^^coZ!) mentioned above. Now, this method 
is not applicable any more. At fourth order, there are 
many eigenvalues, and sometimes they lie close to each 



other and have comparable strength in the response func- 
tion, so that it is not clear which one should be chosen. 

The question arises what is the physical meaning of 
several poles if there is in reality only one damped col- 
lective mode. In order to get a better understanding of 
this question, let us have a look at a simpler example, 
namely a zero-sound wave in a uniform system. For this 
case, comparisons between the moments method up to 
very high order and exact solutions exist in the litera- 
ture [la. Il7j|. In the zero-temperature case, it was found 
|16| that, with increasing order of the moments method, 
the distribution of sharp peaks in the response function 
(i.e., poles just below the real a; axis) converges to the 
continuous spectrum (i.e., a branch cut just below the 
real w axis) of the exact solution of the Vlasov equation. 
Hence, in order to extract the Landau damping from the 
results of the moments method, one has to consider the 
distribution of eigenfrequencies rather than look at their 
imaginary parts. In the case of finite temperature jl7| . 
the collisions provide an additional damping mechanism 
and they lead to complex eigenfrequencies. 

From the preceding discussion it is clear that the fre- 
quency and damping of a mode cannot be obtained from 
the real and imaginary parts of the individual eigenfre- 
quencies given by the moments method, but that one 
has to consider the total response function. This point of 
view is confirmed by the good agreement between the re- 
sponse functions obtained by the fourth-order moments 
method and by numerical simulations in Ref. (9|. 

Besides this theoretical question, there is a more prac- 
tical point one should consider. The idea is that we want 
to compare with experimental data, which were obtained 
by fitting the observed oscillation of the cloud with an ex- 
ponentially damped cosine function. More precisely, in 
the case of the quadrupole mode, the observed oscillation 
is fitted with a function of the form I4I 



Qfitit) = Cie^^ * cos(a;i + ip) + Cse" 



(26) 



while in the case of the scissors mode, the oscillation is 
either fitted with 

Q^7-^(i) = Ce"" cos{ut + ip) (27) 

at low temperature (hydrodynamic regime), or with 



ffit"'-'^{t) 



E^^^ 



-Ffct 



cos{LOkt + ipk) (28) 



fe=i 



at high temperature (coUisionless regime) [l|. So, we will 
determine the frequency and damping rate corresponding 
to our response function Q{uj) by fitting it with Eq. (^5)) 
in the case of the quadrupole mode and with Eq. ((27)) or 
(P5|) in the case of the scissors mode. In the case of a fit 
with two frequencies, we concentrate on the mode with 
the higher frequency. 



TABLE 1: Trap parameters of the Innsbruck experiments. 
Both experiments were done with 600 000 atoms of ®Li in 
the unitary limit {l/kpa = 0) Qj. 



mode 



quadrupole 



iOa:/2TT (Hz) iUy/2TT (Hz) a;z/27r (Hz) 



1600 
1800 



700 
1800 



30 
32 



III. RESULTS 

A. Scissors and quadrupole strength functions 

In Fig. [T] we plot the results for lmQ{uj) obtained at 
second and fourth order for the scissors and quadrupole 
modes at unitarity. The parameters of the trap and 
the number of ^Li atoms are chosen as in Ref. 0], so 
that a comparison with the experimental data is possi- 
ble, see Table HI In the upper panels of Fig. [1] the scis- 
sors response is plotted at various temperatures {T /Tp = 
0.4, 0.6, 0.8). Since mean-field effects are not taken 
into account, the limiting frequencies for the scissors 
mode in the hydrodynamic and coUisionless regimes are 

± w,,, respectively 



^SM = \h^l+^l and ujs^ci± 



W.T 



■Jy, 



[18|. (In the coUisionless regime, two different modes 
can be excited.) In the trap under consideration, these 
frequencies are uis,hd — 1-65 Wr, 0Js,ci- — 0.85 oJr, and 
ijJs.ci+ — 2.17 Wr, where w^ = y/ujx^y is the average radial 
frequency. Let us first analyse the second-order (dashed) 
curves. At T jTp = 0.4, the response is peaked almost 
at u)s,hd'- we are in the hydrodynamic regime. As the 
temperature increases, the peak becomes broader (strong 
damping) and gets shifted towards the higher frequency 
u)s,ci+- At second order, the lower mode at cus.ci- is not 
yet visible at T/Tp = 0.8 since it is still too strongly 
damped. The fourth-order results (full lines) deviate 
more and more from the lowest order ones as the tem- 
perature increases. The most striking feature is that the 
shape itself of the response function is modified by the 
inclusion of the higher-order moments. We also observe 
that at fourth order the lower peak at cl) s.ci- is already 
clearly visible at T/Tp = 0.8. 

In the second row of Fig. [1] we plot the results for 
the quadrupole mode. The limiting frequencies of this 
mode in the hydrodynamic and coUisionless limits are 
^Q,hd = \2.ujr and ijJq^cI = ^ijjr, respectively (again with- 
out mean-field). The second order (dashed) results show 
how the peak moves from the hydrodynamic to the col- 
lisionlcss limit as the temperature increases. Consider 
now the fourth-order (full) lines. At T/Tp = 0.4 and 
T/Tp ~ 0.8 the response shows a clear peak, whose posi- 
tion is however displaced towards higher frequencies, as 
compared to the second-order results. At T/Tp = 0.6, 
the shape of the peak itself is deformed, but again its 
ccntroid is moved towards higher frequencies. This is 
in qualitative agreement with our finding in Ref. (9[ that 
the second-order moments method overestimates the col- 
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FIG. 1: The imaginary part of the scissors (first row) and quadrupole (second row) response function as function of the 
frequency at temperatures T/Tf = 0.4, 0.6, 0.8 (from left to right). The dashed lines represent the second-order results, the 
full ones the fourth-order results. The frequency is in units of the radial trap frequency ojr = ^ujx^y The trap parameters are 
listed in Table H] 



lisional effects, i.e., the second-order result is always too 
close to the hydrodynaniic limit. 



B. Frequencies and damping rates 

In order to make a quantitative comparison of our re- 
sults with the data, we extract from Q{uj) the frequency 
and damping of the mode by fitting the response function 
as explained in Sec. IIIDI The results are shown in Fig. 
[2j Let us first look at the results for the scissors mode 
(first row) . Because of the two different fits at low and 
high temperatures (using one or two damped cosine func- 
tions), there are two curves for each method (second and 
fourth order moments). In the range of T /Tp between 
0.7 and 0.8 we plot both curves in order to show the de- 
pendence on the fit. For T /Tp > 0.8, the gas is closer 
to the collisionlcss regime where two modes are present, 
and we keep only the fit with two damped cosine func- 
tions. For T/Tp < 0.7, we show only the fit with a single 
damped cosine function. Let us now compare the results 
obtained with the second-order (dashed lines) and fourth- 
order (solid lines) moments methods. The most impor- 
tant difference is that the transition from low frequency 
(hydrodynamic regime) to high frequency (collisionlcss 
regime) is shifted to lower temperature by the inclusion 



of fourth-order moments. This was to be expected since, 
as we discussed above, the second-order moments method 
overestimates the coUisional effects. Therefore, for tem- 
peratures below 0.9Tp, the fourth-order frequencies are 
in better agreement with the data than the second-order 
ones. Only at high temperatures, it seems that the 
second-order frequency, which approaches the limiting 
value U!s,ci+ much more slowly, is in better agreement 
with the data. Concerning the damping, there is quite a 
big difference between the second- and fourth-order re- 
sults, but it is not really clear whether the fourth-order 
represents an improvement or not. 

The improvement due to the fourth-order moments is 
more clearly seen in the results for the quadrupole mode 
(second row of Fig. [2]). Again, if the fourth-order mo- 
ments are included, the transition from the hydrody- 
namic to the collisionlcss regime happens at lower tem- 
peratures, which greatly improves the agreement of both 
frequencies and damping rates with the data. But the 
difference between second- and fourth-order calculations 
does not only concern the temperature dependence. This 
can clearly be seen in the right figure, showing the damp- 
ing as function of frequency, so that the temperature 
drops out. In this representation, the curve obtained with 
the fourth-order moments almost passes through the er- 
ror bars of the data, which was by far not the case for 
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FIG. 2: Frequency (left) and damping rate (middle) of the scissors (first row) and radial quadrupole modes (second row) 
as functions of temperature as well as the representation damping vs. frequency (right) which is independent of possible 
uncertainties in the temperature measurement. The points with error bars are the experimental data from Ref. Ql, the dashed 
lines are the second-order results, and the solid lines are fourth-order results. The trap parameters are listed in Table |T1 



the second-order results. 



IV. SUMMARY AND CONCLUSIONS 

In this work, wc determined approximate solutions of 
the linearized Boltzmann equation for collective modes 
of trapped Fermi gases by using the method of phase- 
space moments. Here, we concentrated on the radial 
quadrupole and scissors modes. Contrary to previous 
literature @, S 0, we did not only include the lowest 
(second) order moments which are necessary to describe 
the modes, but also the next (fourth) order. A compari- 
son with a numerical solution of the Boltzmann equation 
[01 showed that the fourth order catches already the most 
important effects missed at second order, e.g., the posi- 
tion dependence of the Fermi-surface deformation. We 
therefore decided to apply this method to realistic cases 
in order to be able to compare with experimental data. 

We showed that, if one includes higher than second- 
order moments, the shape of the response function does 
no longer resemble a single Lorentzian. Therefore, if one 
wants to extract the frequency and damping rate of a 
mode, the result depends on the ansatz for the fit func- 
tion which is used. Our determination of these quantities 
is inspired by the procedure used by the experimentalists. 



In the actual calculation of the moments of the collision 
term, A'^"^^ ^ we used the in-medium cross section as de- 
fined in Ref. [y] . In previous works [j, Q it was found that 
the in-medium enhancement of the cross section strongly 
deteriorates the agreement with the experimental results. 
This conclusion was, however, based on calculations us- 
ing only second-order moments. Since the fourth-order 
contributions reduce the effect of collisions, the effect of 
the enhanced cross section is partially compensated. In 
fact, our new results, including the in-medium cross sec- 
tion and fourth-order moments, are in reasonable agree- 
ment with the data. 

Another application of higher-order moments will be 
to quantify the effects of the anharmonicity of the trap 
potential, including also the mean field. Work in this 
direction is in progress. This might be helpful, e.g., for 
understanding the behavior of the frequencies and damp- 
ing rates at high temperature (note that in Ref. 0] the 
experimental frequencies have been roughly corrected for 
anharmonicity effects by dividing them by the measured 
frequencies of the sloshing mode). For a detailed com- 
parison with the experiment, however, many other effects 
should be accounted for, too. For instance, the measured 
quantity is not the response to a (5 pulse, but the relax- 
ation after the system was adiabatically deformed and 
then suddenly released at i = 0. This results, roughly 



speaking, in an additional factor l/w in the Fourier trans- 
form of the response of the system which can have some 
effect on the fitted frequency and damping rate. In ad- 
dition, the observable measured in the experiment is not 
simply proportional to {xy) , but depends also on the dis- 
tribution in momentum space since the density profile 
in the xy plane is measured after an expansion. With- 
out any doubt, it would be desirable to make a complete 
numerical simulation of the experiment, including the ex- 
pansion phase. 
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and one obtains the useful relation 



dX 
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Using the latter, one can check that 
dX 
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(A5) 

Then, it can easily be seen that the elements of M are 
proportional to {X") eg, n = 2, 4, 6. We choose to express 
them in terms of {x")eq- in trap units, the factors of 
proportionality contain the factor N / P, rational numbers 
and ratios of powers of the trap frequencies. 
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Appendix A: Computation of Aij 

In this appendix we give some details about the com- 
putation of the matrix A defined in Eq. (J12p and on its 
final form. 

As in our practical calculations, we will use trap units, 
i.e., all quantities arc made dimensionless by rescaling 
them by appropriate combinations of the atom mass r/i, 
the average trap frequency uj = {cUx^yUtz)^'^, the har- 
monic oscillator length Iho ~ \/ ^/muj etc. 



The matrix M 

In order to compute Mij defined in Eq. (J13p . it is con- 
venient to define six-dimensional hyperspherical coordi- 
nates. To do this, we must first pass to isotropic spatial 
coordinates, and then to dimensionless ones, so that the 
r and p components can be treated together. We define 

x = lho{uj/i^x)X cos-di 

y = lho{uJ /^y)X s'm'diCos-d2 

z = lyioiGjjujz^X'Bva.-dx sint?2 cos 7^3 
Vx = (l/^;io)-'^sin'i9i sin'i?2 sin'i?3 cosz?4 
Vy ~ (l/^;io)^sin'i9i sin'i92 sin-iJa sin7?4COS(^ 
Vz = (l/^/io)^ sin??! sin'i92 sin-iJa sin'i94 sin(y9 . (Al) 

The volume element becomes d'^rd'^p = X^dXdfl^, and 
its angular part is 

dQ^ = sin* di sin"^ 'i?9 sin^ z?3 sviididi}idi}2d'd'id'didip . 

(A2) 
The integration range is [0,(X)[ for X, [0, 27r] for tp and 
[0,7r] for the di. In these coordinates, the equilibrium 
distribution function reduces to 



feg{X) 



1 



The matrix A*™"" 

Notice that in the case of a harmonic trap, the set 
{(j)i} is closed with respect to the operators in the left- 
hand side of Eq. (j6]), therefore one can find a matrix B 
of coefficients such that 



k=l 



(A6) 



Then, it is clear that the matrix A*™"* defined in Eq. 
(|14p can be written as a matrix product 

where M denotes the matrix calculated in the preceding 
subsection. 

The computation of the matrix B is straight-forward. 
In trap units, its elements are simply given by powers of 
the trap frequencies multiplied by integer numbers. 

The matrix A™" 

To compute the matrix elements Af°", that by defini- 
tion are 



A coll 



^'^^'^■^(r,p)/^/dl7-^'^'P-P^' 



(27r)3 "^'^ '^'7 (27r)3y ~" dn m 

X fejeqlil - /:,)(1 - /:,i)A,„i,[0,] , (A8) 

we follow the method outlined in Refs. [y, [l^. In 
the last equation we have used the compact notation 
^coim = 0(r,p) + 0(r,pi) - 0(r,p') - 0(r,p'i). To 
reduce the number of integrals in Eq. (jASp , we define the 
variables k = p-|-Pi, q = (p — Pi)/2 and q' = (p' — p'j)/2 
(remember that energy and momentum conservation im- 
ply IqI ~ Iq'D- In these variables, one can write 

/eg/egl(l — J en)\^ ~ J eql) ~ 



leqjeqiy^ j eq 
1 



1 



^P{oxy2-p.) _^ I ' 



(A3) 



4 cosh I3{E - n) + cosh /3k • q/2m 
1 
cosh (3{E — fJ.) + cosh /3k • q'/2m 



(A9) 



with E = k'^/im + q^ /2m + V. The factor (j)ilS.\4>j] has 
to be computed and rewritten, as the rest of the inte- 
grand, in these variables, too. Then, we define a rota- 
tion that brings k (identified by the angles 0, ^p) to be 
parallel to the z-axis in momentum space. We define R 
the matrix associated to such a rotation and apply it to 
all momenta: the old coordinates are related to the new 
ones by {px,Py,Pz) = R^^iPa,Pb,Pc), and in particular 
(/Cq, kb, kc) = (0, 0, k). Now the integration over 9, ip can 
be performed analytically, since all the dependence upon 
these variables is in the numerator of the integrand. We 
have thus reduced the number of integrals from eleven 
to nine. Next one defines spherical coordinates for q 
and q': their zenith and azimuth angles are 9c^ ^Pc and 
9'^, ip'^ respectively. Since the dependence upon Lp, ip'^ is 
only in the numerator, we can easily integrate over these 
variables, reducing the integral to a seven-dimensional 
one. Finally, the definition of scaled spatial coordinates 
fi = ^Ti renders the trap potential, and therefore the in- 
tegrand, spherically symmetric in the spatial coordinates: 
the integral is reduced to a five-dimensional one. As a 
result, the elements of A"^"" are proportional, through ra- 
tional numbers and ratios of powers of trap frequencies, 
to terms of the same type of the inverse relaxation time 
1/t defined in Ref. Q. More precisely, now there are 
twelve different terms of this type which are of the form 



1 



• da 



, drP^dkk'^dqq'^^^ I d-fdj'F, 
20TT^mJo ^^ dnj_i ' ' 



1 



cosh I3{E — /i) + cosh I3kqj /2m 

1 
cosh f3{E — /i) + cosh/3fcq7'/2m 



, (AlO) 



i = 1, . . . , 12, where E = k^/8m + q^ /2m, + TOaj^f^/2, 
7 = cos 0c and 7' = cos 6*^. The factors Fi are poly- 
nomials of f^, fc^, g^, 7^, and 7' . In particular, Fi = 
1 + 27^ — 37^7' ^ jSuch that Ji is identical to Is given in 
Eq. (B4) of Ref. [6[ , and is in fact the only non-zero term 
of A"^"" at second order. The coefficients Ji are obtained 
numerically via a Monte Carlo integration and used to 
build A™". 



Appendix B: Calculation of the response function 

In this appendix we describe how the Fourier spectrum 
of a generic observable (x) after a generic perturbation 
V can be obtained. 

First, one has to calculate the vector Ui defined in Eq. 
([TG]) [which is simple in the case V = xy, cf. Eq. ((24|) ]. 
Then, one has to express the expectation value of x in 
terms of the coefficients c^. Supposing that {x)eq = 0, 
the expectation value must be proportional to the c^ and 
one can thus write 



{x)=Y,b^C^ = b^ 



(Bl) 



where we changed to vector notation in the second equal- 
ity, b and c being vectors with components hi and Ci, 
respectively. 

Inverting Eq. pip , one obtains 

(x)(a;) = b^i^iuM + A*™"" + A^^^y^a 

(B2) 

Notice that M, A*™"'^, and yl'^o" are real matrices which 
are independent of w. Now we perform the diagonaliza- 
tion 



M"'(yl'™"" + A""") = PDP- 



(B3) 



with D = diag(Ai, . . . , A„). Since the original matrix 
is real, its eigenvalues X^ are either real or they appear 
as complex conjugate pairs. If we identify the real and 
imaginary parts of the eigenvalues as Xk = Ffc -I- iuik, we 
obtain 



-1 D-l nr-l. 



(x)(w) = b^ P{-iiui + D)-^P-^M 



.^{b'^P)k{P-Hl-^a)u 



cj - cjfc + iFfe 



(B4) 



This reduces to Eq. (|25l) in the special case V = x ^ ^V- 
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